% Matlab code to estimate a flexible mixed logit model in wtp space.
% Written by Kenneth Train, first version Oct 21, 2015, 
% Heather Snell added discounting component April 29, 2019.
 format shortG
 clear all

% Do not change the next line. It sets the global variables.
global NP XMAT PROBS Z NZ NCS NDRAWS NV BETAS ZTYPE CrossCorr YesGPU COEF NROWS discountfun T modelname
% OUTPUT FILE
discountfun = "QuasiHyper";

modelname = 'QH.Thumb';
T=65; %max time periods;
 diary off
 delete QH.Thumb.out
 diary QH.Thumb.out
subset=0; % setting subset = 0 will give full dataset
% subset = 1 has consequential > 1
% subset = 2 has consequential > 3
% subset = 3 has all attributes never ignored
% subset = 4 has all attributes either never ignored or sometimes ignored
filename = 'QH.Thumb.mat';
CrossCorr=1;
NBins=20;
WantHessian=1;

WantBoot=0; %set to 1 to run bootstrapping
NReps=250;
NDRAWS=1000;
NGridPts=6001;

% ZTYPE=1 for polynomials.
%      =2 for step functions
%      =3 for splines
ZTYPE=3;

PolyOrder=8;
NLevels=6;  
NKnots=8;

if discountfun=="Harvey" % 
R_Range=[1.1944 1.2164];
P_Range=[0 12.19413];
WTP_Range=[-9.9824	7.106;
           -7.1937	9.1447;
           -0.8247	1.0361;
           -2.1644	1.9384;
            0.1785	0.2129;
           -0.3655	1.0197;
           -0.0723	0.2805];


elseif discountfun=="HM" % 
    R_Range=[1.4554 1.4678];
    P_Range = [0 11.880];
    WTP_Range=[-9.2509  6.7231;%cap
               -8.7103  10.602;%sw
               -0.5120  0.7104;%buffer
               -1.5578  1.4062;%infra
                0.1930  0.2062;%jobs
               -0.3964  1.0456;%quality
               -1.3590  1.5378];%wlife  
elseif discountfun=="Exponential" % 
    R_Range = [0.490, 0.492];
    P_Range = [0, 9.959];
    WTP_Range=[-10.4355	10.3453;
                -7.2819	10.5761;
                -2.6147	 3.8801;
               -11.4394	 9.7886;
                 0.4169	 0.4461;
               -11.9580	13.9516;
               -16.7970	16.4858];

elseif discountfun=="QuasiHyper" %
    R_Range = [0.4027, 0.4231];
    P_Range = [0.8297, 1.1605]; %actually the beta range 
    WTP_Range=[ 0       13.7575; %price range
               -11.7184	 9.7024;
                -7.7486	 9.5902;
                -2.6359	 3.7321;
                10.0187	 9.4345;
                 0.6754	 0.7174;
                -5.7674	 7.4594;
                -3.0512	 4.2668];
               

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Below should not change %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('discounting.data.mat');  %This loads XMAT with the variables below
FirstRun = 1; %This is used and changed in Bootstapping, but don't change here
Nalts = 3; %number of alts per question
Nquest = 5; %number of questions each person answers
NROWS = size(XMAT, 1);
NP = NROWS/(Nalts*Nquest);
NCS = NP * Nquest; 

% 1. Person number (1-NP)            MUST BE THIS. DO NOT CHANGE.
% 2. Choice situation number (1-NCS) MUST BE THIS. DO NOT CHANGE.
% 3. Chosen alternative (1/0)        MUST BE THIS. DO NOT CHANGE.
% 4. Price 
% 5. CAP
% 6. SW
% 7. Buffer
% 8. Quality
% 9. Jobs
% 10. Infra
% 11. Wlife
% 12. treatment
if discountfun=="QuasiHyper"
    NAMES = {'Rate';  'Beta'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
                'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
else
    NAMES = {'Rate'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
    'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
end

IDPRICE=4;

NV = size(NAMES, 1);
%Do not change the next line
COEF = [R_Range;P_Range;WTP_Range];

ThisSeed=1234;
rng(ThisSeed);

%Do not change the following lines. They calculate the number of Z variables
if ZTYPE==1; 
    NZ=PolyOrder*NV;
    if CrossCorr==1;
        if discountfun=="QuasiHyper"
            NZ=NZ+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
        else 
            NZ=NZ+(NV-2)*(NV-3)/2; %will not correlate price or rate
        end
    end
elseif ZTYPE==2
    NZ=(NLevels-1)*NV;
    if CrossCorr==1;
            if discountfun=="QuasiHyper"
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
elseif ZTYPE==3
    NZ=(NKnots+1)*NV;
    if CrossCorr==1;
            if discountfun=="QuasiHyper"
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
end

%Set the starting values for the coefficients of the Z variables
beta=zeros(NZ,1);  

YesGPU=1;

MAXITERS=5000;

%CONVERENCE TOLERANCE.
% Set =[] to use matlab defaults. default is 1e-6
PARAMTOL=[];
LLTOL=[];
tic;
if discountfun ~= "QuasiHyper"
    CreateDrawsDisWtpProbs;
    CheckRate;
    if check_ok==1
       CreateZ;
       disp(' ');
       disp(['Data setup took ' num2str(toc./60) ' minutes.']);
       disp(' ');
       DoEstimationRate;
        paramhat_orig = paramhat;
        Figure_Output9;
        if WantBoot==1
            Boot
        end
    end
else
    CreateDrawsDisWtpProbsQH;
    CheckRate;
    if check_ok==1
        CreateZQH;
        disp(' ');
        disp(['Data setup took ' num2str(toc./60) ' minutes.']);
        disp(' ');
        DoEstimationRate;
        paramhat_orig = paramhat;
        Figure_Output10;
        if WantBoot==1
           Boot
        end
    end

end

save(filename)
diary off
